function  [D_CS, D_Ext, D_Mis1, D_Mis2, D_Mis, D_SW, m_theta, m_D_E, m_E, sd_theta, sd_D_E, sd_E, corr_theta_D_E, corr_theta_E, D_Mis_SumStats, D_Mis_SumStats_1, D_Mis_SumStats_2, D_Mis_SumStats_3, D_Mis_SumStats_4] = Total_Welfare_decomp_WTP_EEgap_tax_heter(weight,param,param_FE,param_FEdemo,id_ij,J,p_elec,t_pigou,d_ext,elec,lifetime,r,MCPF,MT,price, price_denom, rebate_estar_denom, estar_denom, kwh, prod_denom,id_i, id_j, id_m, id_n, type_id2_3_demo1_denom, AV_demo1_denom,  ice_sc_demo1_denom, type_id2_3_demo2_denom,   AV_demo2_denom,  ice_sc_demo2_denom,type_id2_3_demo3_denom,   AV_demo3_denom,  ice_sc_demo3_denom,type_id2_3_demo4_denom,   AV_demo4_denom,  ice_sc_demo4_denom,type_id2_3_demo5_denom,  AV_demo5_denom,  ice_sc_demo5_denom, type_id2_3_demo6_denom,  AV_demo6_denom,  ice_sc_demo6_denom)
                                                                                                                                                                                                                                                    
    param_FE_k=param_FE;
	param_FEdemo_k=param_FEdemo;

dim_k=length(weight)
%parfor k=1:dim_k
for k=1:dim_k

	param_k=param{k};
    theta_k(k)=param_k(3);

%The social price of electricity is computed as follows: p_elec_ext = elec+p_elec+t_pigou;

%Welfare evaluated at the baseline 
d_ext_J = repmat(d_ext',J,1);
t_pigou_0=t_pigou*0;
[CS_0_i{k}, CS_0_avg(k), E_0_i{k}, E_0_avg(k),  W_0_i{k}, W_0_avg(k), Error_0_i{k}, Error_0_avg(k), ext_0_i{k}, ext_0_avg(k), gvt_0_i{k}, gvt_0_avg(k), SW_0_i{k}, SW_0_avg(k), P_0_decision{k}, P_0_decision_j{k}] = ...
Welfare_decomp_k_WTP_EEgap_parfor(param_k,param_FE_k,param_FEdemo_k,id_ij,J,p_elec,t_pigou_0,d_ext_J,elec,lifetime,r,MCPF,MT,price, price_denom, rebate_estar_denom, estar_denom, kwh, prod_denom,id_i, id_j, id_m, id_n, type_id2_3_demo1_denom, AV_demo1_denom,  ice_sc_demo1_denom, type_id2_3_demo2_denom,   AV_demo2_denom,  ice_sc_demo2_denom,type_id2_3_demo3_denom,   AV_demo3_denom,  ice_sc_demo3_denom,type_id2_3_demo4_denom,   AV_demo4_denom,  ice_sc_demo4_denom,type_id2_3_demo5_denom,  AV_demo5_denom,  ice_sc_demo5_denom, type_id2_3_demo6_denom,  AV_demo6_denom,  ice_sc_demo6_denom);

%Welfare evaluated at the social price of electricity 
[CS_i{k}, CS_avg(k), E_i{k}, E_avg(k),  W_i{k}, W_avg(k), Error_i{k}, Error_avg(k), ext_i{k}, ext_avg(k), gvt_i{k}, gvt_avg(k), SW_i{k}, SW_avg(k), P_decision{k}, P_decision_j{k}] = ...
Welfare_decomp_k_WTP_EEgap_parfor(param_k,param_FE_k,param_FEdemo_k,id_ij,J,p_elec,t_pigou,d_ext_J,elec,lifetime,r,MCPF,MT,price, price_denom, rebate_estar_denom, estar_denom, kwh, prod_denom,id_i, id_j, id_m, id_n, type_id2_3_demo1_denom, AV_demo1_denom,  ice_sc_demo1_denom, type_id2_3_demo2_denom,   AV_demo2_denom,  ice_sc_demo2_denom,type_id2_3_demo3_denom,   AV_demo3_denom,  ice_sc_demo3_denom,type_id2_3_demo4_denom,   AV_demo4_denom,  ice_sc_demo4_denom,type_id2_3_demo5_denom,  AV_demo5_denom,  ice_sc_demo5_denom, type_id2_3_demo6_denom,  AV_demo6_denom,  ice_sc_demo6_denom);

E_0_ik(:,k)    = E_0_i{k};
E_ik(:,k)      = E_i{k};
D_CS_ik(:,k)   = CS_i{k}-CS_0_i{k};
D_Ext_ik(:,k)  = d_ext.*(E_i{k}-E_0_i{k});
D_Mis1_ik(:,k) = (theta_k(k)-1)*elec.*(E_i{k}-E_0_i{k});
D_Mis2_ik(:,k) = theta_k(k)*t_pigou.*E_i{k};
D_Mis_ik(:,k)  = (theta_k(k)-1)*elec.*(E_i{k}-E_0_i{k}) + theta_k(k)*t_pigou.*E_i{k};
D_SW_ik(:,k)   = SW_i{k}-SW_0_i{k};

%We don't consider heterogeneity at the individual (i) level with respect
%to the externality (d_ext), electricity prices (elec), and t_pigou
%emission factors (t_pigou)
%D_CS_avgk(k)   = CS_avg(k)-CS_0_avg(k);
%D_Ext_avgk(k)  = mean(mean(d_ext))*(E_avg(k)-E_0_avg(k));
%D_Mis1_avgk(k) = (theta_k(k)-1)*mean(elec)*(E_avg(k)-E_0_avg(k));
%D_Mis2_avgk(k) = theta_k(k)*t_pigou*E_avg(k);
%D_Mis_avgk(k)  = (theta_k(k)-1)*mean(elec)*(E_avg(k)-E_0_avg(k)) + theta_k(k)*t_pigou*E_avg(k);
%D_SW_avgk(k)   = SW_avg(k)-SW_0_avg(k);

end
    
    E_0_k = sum(E_0_ik,1)/MT;
    E_k   = sum(E_ik,1)/MT;

    D_CS   = sum(weight'.*sum(D_CS_ik,1)/MT);
    D_Ext  = sum(weight'.*sum(D_Ext_ik,1)/MT);
    D_Mis1 = sum(weight'.*sum(D_Mis1_ik,1)/MT);
    D_Mis2 = sum(weight'.*sum(D_Mis2_ik,1)/MT);
    D_Mis  = sum(weight'.*sum(D_Mis_ik,1)/MT);
	D_SW   = sum(weight'.*sum(D_SW_ik,1)/MT);

    %D_CS_avg  = sum(weight'.*D_CS_avgk);
    %D_Ext_avg = sum(weight'.*D_Ext_avgk);
    %D_Mis1_avg = sum(weight'.*D_Mis1_avgk);
    %D_Mis2_avg = sum(weight'.*D_Mis2_avgk);
    %D_Mis_avg = sum(weight'.*D_Mis_avgk);
	%D_SW_avg  = sum(weight'.*D_SW_avgk);

%Compute sufficient stats
m_theta        = sum(weight'.*theta_k); 
m_D_E          = sum(weight'.*(E_k-E_0_k));
m_E            = sum(weight'.*E_k);
sd_theta       = (sum(weight'.*(theta_k-m_theta).^2))^0.5;
sd_D_E         = (sum(weight'.*((E_k-E_0_k)-m_D_E).^2))^0.5;
sd_E           = (sum(weight'.*(E_k-m_E).^2))^0.5;
cov_theta_D_E = sum(weight'.*((theta_k-m_theta).*((E_k-E_0_k)-m_D_E)));
cov_theta_E   = sum(weight'.*((theta_k-m_theta).*(E_k-m_E)));
corr_theta_D_E = cov_theta_D_E/(sd_theta*sd_D_E) ;
corr_theta_E   = cov_theta_E/(sd_theta*sd_E);

%Reconstruct the welfare component based on sufficient statistics
D_Mis_SumStats_1 = sum(elec.*(m_theta-1)*m_D_E,1)/MT;
D_Mis_SumStats_2 = sum(elec.*corr_theta_D_E*sd_D_E*sd_theta,1)/MT;
D_Mis_SumStats_3 = sum(t_pigou*m_theta*m_E,1)/MT;
D_Mis_SumStats_4 = sum(t_pigou*corr_theta_E*sd_theta*sd_E,1)/MT;

%Compare whether the sufficient stats approach matches: D_Mis_SumStats vs D_Mis
D_Mis_SumStats = D_Mis_SumStats_1 + D_Mis_SumStats_2 + D_Mis_SumStats_3 + D_Mis_SumStats_4;

%Compare the sufficient stats approach for the heterogeneity components
%only:
D_Mis_SumStats_2_check = D_Mis1 - D_Mis_SumStats_1;
D_Mis_SumStats_4_check = D_Mis2 - D_Mis_SumStats_3;

end
